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Abstract. To treat the problem of growing protoplanetary disc solids across the meter barrier, we consider a very simplified 
two-component coagulation-fragmentation model that consists of macroscopic boulders and smaller dust grains, the latter 
being the result of catastrophic collisions between the boulders. Boulders in turn increase their radii by sweeping up the dust 
fragments. An analytical solution to the dynamical equations predicts that growth by coagulation-fragmentation can be efficient 
and allow agglomeration of 10-meter-sized objects within the time-scale of the radial drift. These results are supported by 
computer simulations of the motion of boulders and fragments in 3-D time-dependent magnetorotational turbulence. Allowing 
however the fragments to diffuse freely out of the sedimentary layer of boulders reduces the density of both boulders and 
fragments in the mid-plane, and thus also the growth of the boulder radius, drastically. The reason is that the turbulent diffusion 
time-scale is so much shorter than the collisional time-scale that dust fragments leak out of the mid-plane layer before they can 
be swept up by the boulders there. Our conclusion that coagulation-fragmentation is not an efficient way to grow across the 
meter barrier in fully turbulent protoplanetary discs confirms recent results by Brauer, Dullemond, & Henning who solved the 
coagulation equation in a parameterised turbulence model with collisional fragmentation, cratering, radial drift, and a range of 
particle sizes. We find that a relatively small population of boulders in a sedimentary mid-plane layer can populate the entire 
vertical extent of the disc with small grains and that these grains are not first generation dust, but have been through several 
agglomeration-destruction cycles during the simulations. 
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1. Introduction 

The formation of km-sized planetesimals from /art-sized 
dust grains is a long-standing challenge of planet forma- 
tion. The problem is complicated by the interplay of an 
array of different physics - most notably the turbulence 
of protoplanetary discs and the stickin g and collisio n al de- 



struction of solids of d if ferent sizes (ICho kshi et al. 
Weidenschilling & Cuzzil 1 19931 iDominik & Tielena 
Blum 



nschiihng & Cuzzii ID ominik 
& WurmfboOOtlHenning et al.ll2006l) 



1993 



1997 



Young stars are known to receive ma s s from the inner part 



of their circumstellar disc (IBertout et al.L 11988). The cause of 



such accretion is most likely that protoplanetary discs are tur- 
bulent. The source of turbulence must in this connection be 
a Keplerian s hear instability (such as the magnetorotational 
instability, see Balbus & Hawlevi 1 19981) or self-gravity if the 
disc i s massive enough (|Balbus & Papaloizoul 1999 : Gammie , 
20011: lLodato & Ricel |2004 . There are nevertheless signifi- 
cant problems with both these sources of accretion: the ionisa- 
tion fraction of protoplanetary discs at 1-20 AU from the star 
may be too low for the magnetorotational instability to oper- 
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ate dGammie , 1 19961 1 ISemenov et all 120041) . while only a minor 
fraction of discs are expected to be massiv e enough to be grav- 
itationally unstable dBeckwith et al. . 1990h . Thus turbulence is 



often treated as a free parameter in protoplanetary disc mod- 
els, parameterised through an a-value (i.e. turbulent viscosity) 
anywhere from a = 10~ 6 up to as high as a — 0.1. 

Collisions between /itn-sized dust monomers leads to the 
formation of dust agg regates under a range of conditions 
( Blum & Wurml 20001) But larger bodies have poor stick- 
ing p roperties and a lowe r threshold fo r collisional destruc- 



tion dChokshi et all 1 19931: iBenzi |2000|) . The sticking prob- 
lem is especially acute for m-sized boulders. The strength 
of these bodies is very low, while collision speeds, induced 
by the turbulent gas and by the differ e ntial settling a nd ra 



dial drift, are high dVolketall Il980l iMizuno et all Il988 



Weidenschillingl 1997 ). The slightly sub-Keplerian gas acts as 



a constant head wind on the boulders, draining them of angu- 
lar m omentum and causing t hem to drift radially through the 
disc ( Weidenschillingl 1977 ). The drift rate is approximately 
proportional to the radius of the boulders, introducing a differ- 
ential radial drift that can be as strong as 50 m/s in difference 
between bodies of 1 m and 10 cm in size. 
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Modelling the growth of solids in protoplanetary discs re- 
quires solution of the coagulation equation (or Smoluchowski 
equation) that governs the time evolution of a size distribu- 
tion of solids. Pioneering work on the numerical solution of 
the coagulati on equation in a planete simal formation context 
was done by Weidenschillingl ( 1984 ) who found that parti- 
cle growth in turbulent discs is efficient because the relative 
speeds induced by the turbulence give high collision rates. 
Collisional fragmentation nevertheless halts growth when the 
solids reach sizes of a few cm (pebbles). T hese early mod- 
els were improved in Weidenschillingl ( 1997 ) to include many 
more bins in the vertical direction and in the particle radius. 
Considering the formation of comets at 30 AU from the proto- 
Sun, Weidenschilling concluded that coagulation could in prin- 
ciple explain the growth all the way to km-sized planetesimals, 
although for disc models where the turbulence is purely in- 
duced by the sedimentation, so that very high particle densi- 
ties occur in the mid-plane. The sticking efficiency was also 
assumed to be hi gh, in some contrast with later m odels of 
boulder collisions dBenz , l200rj ISchafer et all 12007b . and col- 
lisional fragmentation was kept at a minimum by assuming a 
constant specific kinetic energy threshold for fragmentation. 
In this approach bodies of different sizes may collide at very 
high speeds without d estroying each other (see Appendix F of 
IWeidenschiliing||l997h. 



Dullemond & Dominik (2005) presented simulations sim- 
ilar to Weidenschillind ( 1997 ) and found numerically that the 
size distribution of solids can get into a balance between col- 
lisional fragmentation and coagulation. The first simulations 
to include the full radial extent of a protoplanetary disc were 
presented recently by Brauer, Dullemond, & Henning (2008, 
hereafter BDH). BDH found that the meter barrier is a gen- 
uine problem for planetesimal formation, both because macro- 
scopic bodies destroy each other in catastrophic collisions and 
because radial drift sends macroscopic bodies into the inner 
nebula where they are lost from the planet formation process. 
Radial drift is not only a problem for m-sized boulders: over 
the life-times of protoplan etary discs (millions of years, see 

even mm- and cm-sized pebbles 



Bouwman et al. . 20061) 



e.g. 

have signif icant drift and are emptied from the outer parts of 



the nebula (Takeuch i & Lira 12002; Bra uer et al.L 12007). 



In this paper we isolate the effect of collisional fragmenta- 
tion and subsequent sweep-up of fragments on the growth of 
boulders. We do not solve the full coagulation equation as in 
BDH, but simplify the size distribution to effectively two bins 
- small dust fragments and large boulders - in order to make 
the particle growth tractable in a 3-D simulation with magne- 
tised time-dependent turbulence. 

The paper is structured as follows. In $2] we describe our 
simplified two-component model of boulders and fragments in 
detail and find an equilibrium solution to the dynamical equa- 
tions. In O we describe the numerical simulations that will be 
used to evolve the dynamical equations. A local corotating box 
is considered, and turbulence is produced by the magnetorota- 
tional instability. Dust fragments are treated as a passive scalar, 
while the boulders are treated as individual superparticles with 
two internal degrees of freedom: the number density of actual 
particles inside each superparticle and the average radius of the 



constituent boulders. In the next section, §H] we present simu- 
lations where the dust fragments are not allowed to leave the 
boulder layer. The system quickly evolves towards the equilib- 
rium state found in §|2] with rapid growth of the boulders of a 
few mm per orbit. In §|5]we briefly turn to the analytical model 
again and see how the diffusion of dust fragments from the 
boulder layer affects the equilibrium solution. The turbulent gas 
transports fragments out of the mid-plane, and thus the growth 
rate of the boulders is reduced. This is confirmed in computer 
simulations where boulders lie in a thin layer around the mid- 
plane, presented in §|6] The dust fragments spread quickly out 
of the mid-plane and the growth rate of the boulders is reduced 
by a factor of 10, so that the growth can no longer compete with 
the radial drift. The whole disc is filled with dust fragments 
that are formed in the thin mid-plane layer. We conclude on 
our results in §0 and speculate about ways to get coagulation- 
fragmentation growth back on track. We find that radial drift, 
not collisional fragmentation, is the more serious problem for 
coagulation-fragmentation growth and discuss processes to re- 
duce or stop radial drift in actual protoplanetary discs. 

2. Coagulation-fragmentation model 

We consider a simple two-species model of solids in a proto- 
planetary disc. Species 1 consists of tiny dust grains with mass 
mi and number density n\, while species 2 consists of macro- 
scopic boulders with mass ra 2 and number density n 2 - Here 
"dust grains" are defined as being so small that they couple 
to the gas on a time-scale that is much shorter than an orbital 
time, while "boulders" are solid bodies with sizes somewhere 
between 10 cm and 10 m. 

We assume for simplicity that 

1. Collisions between the tiny grains are insignificant com- 
pared to the sweep-up of the grains by the boulders. 

2. Collisions between a boulder and a dust grain always lead 
to the incorporation of the grain into the boulder. 

3. Collisions between the boulders lead to a complete destruc- 
tion of the colliding bodies. The entire mass then ends up 
in tiny grains. 

Thus we do not treat the problem of how to form boulders in 
the first place, but focus on how they grow by sweeping up 
dust. We also ignore effects like cratering in our treatment of 
boulder-dust collisions. Although the actual growth of boul- 
ders across the metre-barrier will likely involve a combination 
of many different aspects of collision physics (and also self- 
gravity), we will in this paper instead aim at gaining insight 
into the pure problem of collisional fragmentation of equal- 
sized boulders and sweep-up of small dust grains in a turbu- 
lent environment. Assuming an impact strength of zero for the 
boulders and perfect sticking between boulders and dust grains 
may not be entirely realistic, but this allows us to simplify our 
model greatly. 

We can write up the dynamical equations for the number 
densities of dust and boulders and for the mass of the individual 
boulders, 



dni 2 m 2 

— = - nin2(TnVn + H 2 cr 22 V22 — 

at mi 



(1) 
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dn 2 2 

— = -« 2 cr 2 2V 2 2 , 

dm 2 



dt 



0"l2Vi 2 pi 



(2) 
(3) 



Here <x 12 and o"2 2 are the collisional cross sections for boulder- 
grain and boulder-boulder collisions, respectively, while v\ 2 
and V22 are the corresponding collision speeds. We furthermore 
introduced the bulk density of dust fragments p\ = n\ni\. We 
assume next that the grains and the boulders are spheres with 
radius ci\ and a 2 , respectively, and that ai <sc a 2 . The dynamical 
equation for m 2 can then be turned into a dynamical equation 
for the radius a 2 , 



da 2 
dt 



Pi 
4p. 



V12 : 



(4) 



with p. referring to the material density of the solids. 

2A. Equilibrium limit 

Eq. (HJ has the equilibrium solution 

^=4^, (5) 

P2 V12 

where we define the bulk densitie^] pi = n\m\, p 2 = n 2 m 2 
and set cr 22 ~ 4<x 12 under the assumption that the contribution 
of the small grains to the cross section cr\ 2 is vanishing. We 
have also assumed, by setting the collisional radius of a boul- 
der to twice its physical radius, that all impact parameters lead 
to destruction, even if the boulders collide at a small grazing 
angle. Thus if V\ 2 ~ v 22 the system tends towards an equilib- 
rium where the small grains in total contain four times more 
mass than the boulders. We show in Appendix [A] that any per- 
turbation to the equilibrium solution will decay on a collisional 
time-scale, so that Eq. © constitutes a (both linearly and non- 
linearly) stable solution to the coagulation-fragmentation prob- 
lem [Eq. ©]■ 

Inserting Eq. (O into Eq. (0]i yields the evolution of the 
boulder radius in the equilibrium state as 



da 2 
~dt 



Pl+2 



V'22 



p. 1 + 4V22/V12 



(6) 



Here we have introduced the bulk density of solids pi + 2 = 
pi + p2 which is constant in time in absence of evaporation and 
condensation processes. We show in Fig.[TJthe dependence of 
the radius growth on the collision speeds V\ 2 and V22- It is clear 
from Eq. © that the radius growth depends only on either v\ 2 
or V22 in the two limits of \ 22 jv\ 2 (i.e. zero and infinity). The 
dividing line at v 22 /vt 2 = 1/4 is indicated with a black line in 
Fig.ffl 

Eq. (O implies that there is a linear growth of the boulder 
radius with time, although only under the assumption that pi+2, 
V12 and V22 are independent of particle size (we shall include 
the full dependence of these parameters on the particle size in 



1 We use the term "bulk density" throughout this paper to refer to 
the total mass of solid material in a given volume divided by the vol- 
ume, i.e. including the void between the solids. 



0.20 



0.15 



0.10 



0.05 



0.00 
0.00 




0.05 



0.10 



0.15 



0.20 



Fig. 1. Contour lines of the radius growth of boulders, a 2 , 
as a function of the collision speed between boulders and 
fragments, vn, and between boulders and boulders, V22- Two 
regimes are divided by the black line: for V22/V12 <K 1/4 the 
radius growth depends only on V22, whereas the radius growth 
depends only on v\ 2 in the limit V22/V12 » 1/4. Note that the 
normalisation of collision speeds with sound speed is an arbi- 
trary choice. 

the numerical simulations presented in §0 and §|6j. For typical 
values ofpi + 2/p. = 10 -11 , V12 = 25 ms -1 and V22 = 10ms _1 , 
relevant in a sedimentary mid-plane layer of solids at r — 5 AU 
of a moderately turbulent minimum mass solar nebula model 
with a turbulent viscosity of a — 10~ 3 (see ^3] for a definition 
of a), the growth rate is a 2 — 4 x 10 _I 1 m s , or 1 .2 millimeters 
per year. Around 2,000 years are then needed to grow from 30 
cm, the size for which radial drift is the fastest, to 3 m in radius, 
which is so loosely coupled to the gas that radial drift is no 
longer a problem. In the absence of collisional fragmentation, 
on the other hand, the sweep-up will end after the boulders have 
incorporated all the small grains. Considering a fixed number 
density of boulders n 2 , we can write the ratio of the particle 
radii a 2 and a' 2 for two different mass densities p2 and p' 2 as 



pi 



(7) 



Setting p' 2 = pi + p2, it is seen from Eq. (O that if 4/5 of the 
dust mass is originally in small grains, then the sweep-up of 
those grains by the boulders can only lead to a moderate in- 
crease of approximately 70% in the average boulder radius. 
Only when collisional fragmentation is included can the boul- 
ders grow larger than that, because the reservoir of grains to 
sweep up will be continuously replenished. 

For the case v\ 2 ~ v 22 , Eq. © simplifies down to 



da 2 



1+2 



5 p. 



-v 22 . 



(8) 
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One can consider yet another special case of Eq. © where 
the collisions between boulders and grains happen at a much 
higher speed than the collisions between boulders and boul- 
ders, V12 » V22- This is relevant if the boulders migrate ra- 
dially inwards due to a radial pressure gradient in the gas 
(IWeidenschillingll 1977b . The drift speed can approach 10% of 
the sound speed, easily an order of magnitude higher than the 
turbulent gas motions that cause the collisions between boul- 
ders. Thus sweep-up works much more efficiently than frag- 
mentation, and Eq. (O changes to 



~dt 



Pl+2 



■V22 ■ 



(9) 



which is five times faster than Eq. ([8j. The approximation that 
V22 ^ V12 may nevertheless be unachievable, even if turbulent 
motion is weak, since shape effects will in duce differen tial ra- 
dial drift even between equal-mass bodies dBenzL l2000h . 




10 20 30 40 50 



2.2. Timescale to reach equilibrium 

Starting from a state where an equal amount of mass is present 
in dust and in boulders, the coagulation-fragmentation equilib- 
rium is reached when a significant fraction of the boulders have 
undergone collisions (and the following fragmentation). From 
Eq. (f2]i the time it takes to get to equilibrium r eq is given by 



1 



«20"22V22 ■ 



(10) 



Assuming that the boulders are spherical, the number density 
«2 and the collisional cross section o"22 can be written in terms 
of the solid radius a%, yielding 



p.ai 

3p2V22 



(11) 



When a.2 < (9/ 4) A, where A is the mean free path of the gas, 
the friction force is in the Epstein regime (see Appendix El for 
a discussion of the validity of the Epstein regime for the boul- 
ders). Here the friction time can be written as 



Tf = 



a 2 p. 

CsPg 



(12) 



with p g denoting the gas density and c s the sound speed. 
Inserting this expression for the friction time in Eq. ( fTTT i yields 



3e 2 V22/c s 



(13) 



Here Qk = is the Keplerian angular frequency of the 

disc, a measure of the local dynamical time-scale at a given ra- 
dial location in the disc, and ex is the ratio of the bulk densities 
of boulders and gas. We define a dimensionless friction time 
through the Stokes number St as 



St = Q K T f . 



(14) 



Thus the time it takes to reach a coagulation-fragmentation 
equilibrium is independent of the actual density of the sur- 
rounding gas. For a given Stokes number, any decrease in the 



Fig. 2. The time evolution of the mass density in dust grains, p\ , 
relative to that of boulders, p2- We have assumed an initial dust- 
to-gas ratio of 0.3 for both dust and boulders, collision speeds 
of V12 = V22 = 0.02c s and a solid density of p./p g = 10". The 
timescale f eq for approaching the equilibrium value p\ /p 2 = 4 
(dashed line) is indicated with vertical lines. 



gas density must be balanced by a similar decrease in the ra- 
dius of the boulders, leading to an increase in the number den- 
sity that balances out the decrease in collisional cross section in 
Eq. dTOb . The presence of the sound speed in Eq. fl3i also does 
not affect the timescale, since the turbulent collision speed v'22 
must scale with the sound speed as welS 

The time evolution of the mass ratio pi/p2 between dust 
grains and boulders is shown in Fig. [2] We have integrated the 
0-D coagulation-fragmentation equations [Eqs. (Q~|), (f2]i and (|4|i] 
with an initial dust-to-gas ratio of 0.3 for both species (relevant 
in a mid-plane layer of solids in equilibrium between sedimen- 
tation of turbulent diffusion), collision speeds of V12 = V22 = 
0.02c s and a solid density of p./p g = 10 11 . The approach to 
the equilibrium state P1/P2 = 4 happens on the equilibrium 
timescale r eq , given by Eq. ( fTTT i. For St = 0.5 the equilibrium 
timescale is around five orbits, increasing proportional to the 
Stokes number. The onset of fragmentation should depend on 
the radius of the boulders rather than on the Stokes number, 
so the critical Stokes number for which collisional fragmenta- 
tion gets important depends on the radial location in the disc. 
Around the location of Jupiter in a minimum mass nebula a 
St = 1 particle has a radius of approximately 30 cm (see 33. 4K 
for which collisional fragmentation should already be substan- 
tial. Since the coagulation-fragmentation equilibrium is stable 
(see AppendixlAl. the state will stay at the equilibrium once it 
is reached. 

It is not strictly necessary to be near the equilibrium value 
of P1/P2 to have efficient growth of the boulders by sweep-up. 



2 This is strictly not the case in the presence of magnetic fields, 
where the local Alfven speed gives a second velocity scale, but we 
shall ignore that complication here. 
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Run Resolution Particles Leaking p./p g ^il^g ^il^s a ilH a Simulation time 



A 


64 3 


2.0 x 10 s 


No 


10" 


0.3 


0.3 


10-" 


io- 3 


20or orb 


B 


64 3 


2.0 x 10 6 


No 


10 11 


0.3 


0.3 


10-" 


io- 2 


20or orb 


C 


64 3 


2.5 x 10 5 


Yes 


10" 


0.01 


0.01 


10 !1 


io- 3 


40or orb 


D 


64 3 


2.5 x 10 5 


Yes 


10 11 


0.01 


0.01 


io-" 


10- 2 


40or orb 


E 


128 3 


2.0 x 10 6 


Yes 


10 11 


0.01 


0.01 


io-" 


io- 3 


20or orb 


F 


128 3 


2.0 x 10 6 


Yes 


10" 


0.01 


0.01 


io-" 


io- 2 


20or orb 



Table 1. Simulation parameters. The box size is fixed at (1.32#) 3 in all runs. The initial column densities E\ and E2 are set in 
runs A and B to mimic the density in a sedimentary mid-plane layer; in the simulations with vertical gravity on the boulders 
(runs C-F), where dust fragments can leak freely out of the boulder layer, we set the column densities to a more canonical value 
of 0.01 for each component. 



The radius of the boulders grows proportionally to p\, accord- 
ing to Eq. |@), so just maintaining the reservoir of grains is 
already an achievement of the collisional fragmentation. As p\ 
increases with time, the boulders will grow faster and faster 
until finally reaching the growth speed given by Eq. (|6]). 



frequency of lShakura & Sunvaevl (11973b . v t 
obtain the a-value of the disc through 



a — (3nY 



ac 2 Q K , we 



(16) 



3. Simulation set up 

Next we will solve the coagulation-fragmentation equations 
numerically in a three-dimensional time-dependent turbulent 
flow. In this section we describe the dynamical equations and 
the adopted protoplanetary disc model. We let turbulence arise 
throu gh the magnetorotational instability (Ba lbus & Hawlev , 
19911) wh ich operates when the gas is sufficiently ionised 



( lGarnmielll996tlSemenov et al.L 120041) . Typically the most un 



stable wavelength of the magnetorotational instability is around 
one gas scale height, with a subsequent energy cascade to 
smaller scales that approximately obeys a Kolmogorov-law 



(Ha wley et all 119951) . The saturated state of the magnetorota- 



tional instability (which we will refer to as magnetorotational 
turbulence) is characterised by an outwards transport of angular 
momentum through positive Reynolds and Maxwell stresses. 
In shearing box simulations the measured a-value ranges from 
10~ 3 (with zero net flux fiel d) to above 0. 1 (for/3 = P g!lfi /P mas , = 



400 net vertical field, see lHawlev et al.L 11995b . Numerically, 



magnetorotational turbulence has the great advantage that it is 
relatively easy to produce and sustai n in local box simulations 



for hundreds of disc rotation periods (iBrandenburg et al.L 1 1995 



Hawlev et al.. 1995). 



Gullbring et al 



Sicilia-Aguilar et al 



(11998b and more recently 
(120041) measured the accretion lumi- 
nosities of T Tauri stars and translated the measurements 
into mass accretion rate M. Typical estimated values of the 
mass accretion rate lie in the interval M = lO _9 ' _7 M yr _1 . 
Coupling the mass accretion rate with a disc model y ields the 
turbu lent viscosity of the disc, v t , through the relation ( IPringle , 
1981b 



(15) 



where 27 is the column density of gas and solids. Making use 
of the non-dimensionalisation with sound speed c s and angular 



For the minimum mass solar nebula a = 10~ 4 ' "~ 2 from typical 
mass accretion rates. The turbulent viscosity can be approxi- 
mated as 



i 2 

•*rms 1 



(17) 



where r e dd y is the eddy turn over time and « rms is the turbulent 
rms speed. Assuming r e[ jdy ~ £? K ' , du e to the dominating effec t 



of the Coriolis force at large scales ( Weidenschillind, 1984 ) 
one obtains a — (u ms /c s ) 2 . However one must be careful when 
translating a into M lms this way, since a normally refers to the 
turbulence's ability to diffuse the main Keplerian differential 
rotation, and instabilities that are not Keplerian shear instabili- 
ties are often associated with a negat ive g-value (such as con- 
vection or streaming insta bility, see Ryu & Goodman! 1992 : 
lYoudin & Goodmanll2005b . For magnetorotational turbulence, 
Eq. (T7\ nevertheless holds relatively well (see Table©. 

In this paper we shall focus on two values for the viscosity 
which we believe are most relevant (based on the observed ac- 
cretion rates): low viscosity with a = 10~ 3 (arising in zero net 
flux simulations) and high viscosity with a = 10~ 2 (the result 
of simulations with a w eak @ = 20000 vertic al magnetic field). 
We use the Pencil Code ( Brandenburg! 2003 ) to solve the equa- 
tions of ideal magneto hydrodynamics, as described in detail in 
Johansen et al.l (12006b . For simplicity we ignore vertical strat- 
ification of the gas and model a local corotating shearing box 
with side lengths 1.32//, where H = c s /£2k is the scale height 
of the gas. Our coordinate frame is oriented in such a way that 
the A-axis points outwards along the radial direction, the >'-axis 
points along the main Keplerian flow, while the z-axis points 
vertically out of the disc in the direction of the Keplerian fre- 
quency vector £2 K . The simulation parameters are written in 
Table Q] The initial condition for the solids is explained in 33.41 

3.1. Sweep-up 

Boulders are treated as individual particles, each with a unique 
position and velocity vector. The boulders feel a drag force 
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from the gas, described in detail in Appendix [B] but for sim- 
plicity we assume that gas feels no drag from the boulder com- 
ponent. Each particle represents a huge number of actual boul- 
ders, hence we refer to them as superparticles. The dust com- 
ponent is treated as a passive scalar: the velocity field is set 
equal to that of the gas, so that only a continuity equation must 
be solved, but with additional source and sink terms (due to 
destruction of boulders and sweep-up) as described below. 

The superparticles are given two internal degrees of free- 
dom - the number density of actual boulders inside each parti- 
cle hi and the average radius a, of the constituent boulders. A 
superparticle is allowed to change the radius of its boulders, a,-, 
by sweeping up dust grains. The dynamical equation for a, is 



ddi _ £\p g 
~dt ~Jp~. 



(18) 



Here e\ is the dust-to-gas ratio of the dust grains, p g is the den- 
sity of the gas, p. is the material density of the solids, and Vo; is 
the relative speed between superparticle i and the gas in its grid 
cell (the dust grains are so coupled to the gas that this is the 
same as the collision speed between boulders and grains). The 
density of the dust grains is depleted at the same time according 
to the evolution equation 



dei 

dt 



(19) 



with hi denoting the number density of boulders in the super- 
particle i. This number does not change in a sweep-up process. 

3.2. Collisions between boulders 

We identify all superparticles that reside in the same grid cell 
as colliding. Collisions between boulders in the superparticles 
; and j happen at the rate 



Cij = O-ijtlitljVij , 



(20) 



where <Tjj is the collisional cross section and vy is the collision 
speed. We assume spherical particles with crij = 7r(a,- + aj) 2 . 

Boulder collisions are assumed to always lead to total de- 
struction of the colliding bodies. The internal number densities 
of the colliding superparticles i and j change as 



dhj 
~dt 

dt 



= -en, 



(21) 
(22) 



which has to be considered for all combinations of ; and j in 
each grid cell. The total mass that is lost in destructive colli- 
sions is subsequently transferred to the dust component. Here 
the dust-to-gas ratio increases as 



dt\ 4 3 3 



One can think of many improvements to these simplified 
coagulation-fragmentation evolution equations, but we believe 
that it is enlightening to consider the most simple dynami- 
cal equation system that displays coagulation-fragmentation 



growth. We shall also compare our results to the advanced mod- 
els presented in BDH where a size distribution of solids is con- 
sidered and where the fragmentation model is much more ad- 
vanced and show that our results are in relatively good agree- 
ment with this more advanced model. 

One may suspect that the sedimentary mid-plane layer 
could be dense enough to have a gravitational influence on the 
produced fragments. The gravitational acceleration of a homo- 
geneous, infinitely extended mid-plane layer with density pro- 
file p p (z) = [2p/( V2^// p )]exp[-z 2 /(2i/ p 2 )] is 



-2nGZ v exf 



V5ff, 



(24) 



p i 



This acceleration is 3-4 orders of magnitude smaller than the 
vertical gravity from the central star (g z = -Q^z). Thus the 
self-gravity of the sedimentary mid-plane layer can have no in- 
fluence on the escape rate of fragments and we shall ignore the 
effect of self-gravity in this paper. 

3.3. Radial drift 

The accretion process causes the gas pressure in protoplanetary 
discs to fall with radial distance from the young star. We can 
write the global radial pressure gradient acceleration as 



1 dP 

p g dr 



c 2 d\nP 
r d\nr 



HdlnP 



dlni 



(25) 



Here H/r is the disc aspect ratio and c s is the sound speed. 
The balance between Coriolis force and global pressure gradi- 
ent gives the gas orbital velocity relative to the Keplerian mo- 
tion as 



LI 



(subK) 

y 



IHdlnP 
2 r <91nr 



(26) 



It is common to define the pressure gradient parameter r\ as 
( INakagawa et al. . 119861) 



d\nP 



8 In 



(subK) 



(27) 

t/vk- Here vk = QyJ is the Keplerian orbital 



giving u 
speed. 

The boulders do not feel the global pressure gradient and 
would orbit with the local Keplerian speed in absence of gas. 
The head wind of the slower moving gas drains the boulders 
of angular momentum and imposes a flux of boulders towards 
smaller r. The equilibrium radial drift velocity of the boulders 
is given by dWeidenschillina,ll977l:lYoudin & Johansenl 120071) 



v x 



2tjv k 



Q k t { + (Q k t { )- 



(28) 



Typical values o_f_ wvk lie between 0.02c s and 0.1c s 

1986). In this paper we assume that v x = 



(23) (INakagawa et al 



- 0.05c s for marg i nally coupled boulders with £? K Tf = 1. As 



Johansen et al.l ( 120061) we apply the global pressure gradient 



force directly on the boulders instead of on the gas. This sim- 
plified treatment of radial drift is valid as long as the drag force 
from the boulders on the gas is ignored. 
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3.4. Units and initial condition 

We adopt a dimensionless unit system in our corotating box by 
setting the sound speed c s = 1, Keplerian frequency = 1 
and mid-plane gas density p s (z = 0) = 1 (for computa- 
tional simplicity we ignore gas stratification, so the gas den- 
sity is approximately one everywhere in the box). Thus gas and 
particle velocities are measured in units of the sound speed, 
while length and particle radius is in units of gas scale heights 
H = c s /I?k = 1 ■ F° r the magnetic fields we set the vacuum per- 
meability /io = 1 ■ We stress that the role of magnetic fields in 
our model is to tap into the Keplerian motion and release gravi- 
tational energy as turbulent kinetic energy, which is transferred 
from the gas to the solids by drag forces, but that dust grains 
and boulders are otherwise unaffected by magnetic fields. 

Physical parameters for the different runs are written in 
Table [TJ The initial solids-to-gas ratio of boulders and frag- 
ments is set to 0.3 for both species in the two runs where 
dust fragments are not allowed to diffuse out of the boulder 
layer (runs A-B), to mimic the density in a sedimentary mid- 
plane layei0 , whereas the ratio between solids and gas col- 
umn densities is set to the more canonical 0.01 in the simu- 
lations with vertical gravity acting on the boulders (runs C-F). 
The dynamical equations of coagulation-fragmentation further- 
more depend on internal properties of the solids: the radius of 
the boulders a 2 and the material density of the solids p.. These 
must be defined in code units. We initially set a 2 = 10~ n and 
p, = 10 11 , giving an initial Stokes number of unity through 
St 2 = Q k t 2 = (fla/fl)(p./Pg) = 1. 

At r = 5 AU in the minimum mass solar nebula with 
sound speed c s = 5 x 10 4 cms _1 , Keplerian frequency = 
1.7 x 10~ 8 s~' and gas column density £ g = 150gcirT 2 our 
dimensionless model corresponds to a gas scale height of 
H = 3 X 10 12 cm, a mid-plane gas density of p g (z = 0) = 
2 x 10~ n gem -3 , a material density of p. = 2gctrT 3 and an 
initial boulder radius of a 2 = 30 cm. 



4. Results neglecting leaking 

We first treat models where the dust fragments are not allowed 
to leave the boulder layer (runs A and B in Table [TJ to test the 
validity of the analytical model described in §f2] We mimic the 
physical conditions in a sedimentary mid-plane layer of boul- 
ders by setting the densities of boulders and dust grains artifi- 
cially high (both components are given a solids-to-gas ratio of 
0.3). In the following sections, §|5]-[U we treat tne more realistic 
case where boulders lie in a thin layer around the mid-plane and 
where the dust fragments can leave this layer due to turbulent 
diffusion. We caution the reader already now that allowing the 
dust fragments to escape from the sedimentary mid-plane layer 
will make prospects to cross the meter barrier by coagulation- 
fragmentation much more negative than they appear in this sec- 
tion. 



In Fig. [3] we show snapshots of the density of fragments 
(relative to the local gas density) as a function of time. The tur- 
bulence has been given 20 orbits to develop before the sweep- 
up and fragmentation terms are turned on, to avoid the initial 
condition having any influence on the results. After one orbit 
(at t = 2ir 01 b) fragments have formed in collisions between 
boulders. The fragments are continuously mixed by the turbu- 
lence, and after a few orbits the disc reaches a state where the 
fragments are very well-mixed with the gas. This state is pre- 
ferred even though the boulder component is not homogeneous, 
because the collision time-scale of boulders is much longer 
than the diffusion time-scale. Considering an overdense region 
of size A, the time-scale for fragments to diffuse out of this re- 
gion is fdiff = A 2 /D t , where D t is the diffusion coefficient, while 
the collisional time-scale is f co n = ^/(jc^cr^v^). The ratio of 
the two time-scales is 

3(A///) 2 /<5 



fdiff 



(29) 



fcoll £ K Tf(v 22 /c s ) l {pilp s ) 1 

Here we have used the parametrisation D t = 5c 2 Q^ and the 
Epstein friction time Tf = a 2 p./(c s p g ). Using typical values for 
run B, A/H = 0.1, 6 = 10~ 2 , Q K T f = 1, v 22 /c s = 0.07 and 
Pi/Pg = 1 in the overdense regions, yields fdiff/fcoii ~ 0.2. Thus 
the fragments have plenty of time to escape the overdense re- 
gions before they are swept up by the boulders there, leading to 
an almost homogeneous spatial distribution of dust fragments 
at t = 25T or b in Fig. [3] 

In Fig.@]we show the evolution of the friction time of boul- 
ders (first panel), the relative speeds between fragments and 
boulders vi 2 and between boulders and boulders v 22 (second 
panel), the ratio of total fragment mass to total boulder mass 
M\ IM-2 (third panel), and the analytical expectation for M\ /M 2 
(fourth panel) based on Eq. (|5). The growth of boulders by 
sweep-up is very efficient and allows growth to Q^ji > 10 in 
100-200 orbits, around the same time-scale as the radial drift. 
At this size (approximately 10 meters) radial drift is insignif- 
icant and the boulders are no longer in any risk of being lost 
to the inner part of the disc. The variation in particle radius at 
a given time (not shown in Fig. |4j) is generally within 10% of 
the average value, due to the strong coupling between different 
regions of the flow by turbulent diffusion. 

Average values of collision speeds and boulder growth rates 
for all the simulations are written in Table [2] The growth rate 
is almost twice as high in run B (with a = 10~ 2 ) than in run 
A (a — 10~ 3 ) because of the higher collision speeds in the 
strongly turbulent case. 

4.1. Clumping 

There is a discrepancy of around a factor 2-3 between the an- 
alytical expectation of Mi /M 2 in Fig. |4](panel 4) and the mea- 
sured value (panel 3). This is likely due to the fact that the boul- 
der layer is not homogeneous because the pa rticles are con 



centra ted in high pressure regions of the gas ( IJohansen et al 



3 Actually it is not very realistic to have such a high density of small 
dust grains in a sedimentary mid-plane layer to begin with. Therefore 
we also ran variations of runs A-B with all the mass initially in the 
boulder component, but found essentially the same results. 



2006). Fragments are primarily produced in the overdense re- 
gions, but they quickly mix in with the gas (which is approx- 
imately isodense since the turbulence is subsonic). Thus the 
high density regions must produce more fragments to keep up 
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Fig. 3. Time series of the concentration of small dust grains for the run with a = 0.01 and no vertical gravity on the boulder 
component (run B). The boxes are oriented with the radial x-axis to the right and slightly up, the azimuthal y-axis to the left and 
up, and the vertical z-axis directly up. The time t is given in local orbital periods - coagulation and collisional fragmentation are 
turned on after 20 orbits to avoid any effect of the initial condition on the evolution of the solids. Fragments immediately form in 
catastrophic collisions between boulders. Even though the boulder component is clumpy, evident from the inhomogeneous initial 
production of fragments, the dust fragments eventually spread out evenly through the box. This homogenisation is an effect of 
the turbulent diffusion time-scale being much shorter than the collisional time-scale. 



with the turbulent diffusion, and this increases the overall ratio 
of fragments to boulders. 

One could have thought that local overdensities in the boul- 
der layer would allow for enhanced collisional fragmentation 
and thus faster growth out of the radial drift regime. This ap- 
pealing picture nevertheless turns out to be incorrect. To see 
the effect of clumping we imagine collecting the material from 
N grid points into one single grid point, or equivalently to take 
the boulders from a volume V and press them together in the 
volume V/N. The total mass density is N(pi+2)- The single grid 
point that contains boulders must fulfil Eq. © in order to be in 
equilibrium. This leads to the equation system 



Np\ +p 2 = N{pi +2 ), 
Pi-ACpi = 0, 



(30) 
(31) 



where we define f = V22M2 an d assume that the diffusion time- 
scale is much shorter than the collisional time-scale so that p\ 
will be constant among all grid cells. The solution to the alge- 
braic equation system is 

N 

Pi = * r , 1 1, a ~ <Pi+2>> (32) 



N +1/(40 
N 



P2 



-<Pl+2> ■ 



(33) 



For N = 1, corresponding to no clumping, we recover the 
usual pi = (4/5)<pi +2 > and p 2 = (l/5)<pi +2 > for f = 1. 
For N — » 00 the expressions tend towards p\ = (p\+z) and 
Pi — [l/(4£)](pi+2). Thus clumping has little or no effect on 
the equilibrium density of dust fragments, which is the cru- 
cial parameter that determines radius growth. The overdense 



Johansen et al.: The growth and destruction of preplanetesimals 



9 




t/T mb t/T mb 

Fig. 4. Growth and fragmentation of boulders as a function of time t, measured in orbits, for two different strengths of the 
turbulence. The first panel shows the friction time of the boulders: growth from the initial Q-^Tf — 1 to beyond Q^Tf = 10 
occurs readily for both weak and strong turbulence, although the stronger turbulence helps growth by producing more fragments 
(note that the Epstein drag law is formally not valid beyond Q^ts ~ 7, see discussion in Appendix [B] but since this state can 
already be considered as having crossed the meter barrier, we have ignored the complications of switching to a Stokes drag law). 
The second panel shows the relative speeds between boulders and fragments (V12) and between boulders and boulders (yyi). The 
boulder collision speed decreases with time as the particles grow and decouple from the gas, whereas the relative speed between 
grains and boulders is set primarily by the turbulent motion of the dust grains and thus stays approximately constant. The total 
mass of the fragments, M\, reaches 2-3 times the mass of boulders, M2, in the equilibrium state (panel 3), somewhat higher than 
the analytical expectation (panel 4) based on Eq. (|5}, but that is likely due to the fact that the boulder density is not homogeneous: 
fragments are created where the boulder density is high, but quickly escape these regions by turbulent diffusion, leading to an 
increase in the total amount of dust. 



Run Res 


Leaking 


a 




V12 


V>22 




a 2 


frecyc 


A 64 3 


No 


10 


3 


0.067 


0.023 


2.82 


0.057 


3.1 


B 64 3 


No 


10" 




0.12 


0.071 


2.54 


0.089 


2.5 


C 64 3 


Yes 


10" 


3 


0.066 


0.014 


94.0' 


0.0036 


803.5 


D 64 3 


Yes 


10- 


2 


0.12 


0.054 


57.6 


0.0064 


306.4 


E 128 3 


Yes 


10- 


3 


0.063 


0.010 


58.0* 


0.0029 


640.5 


F 128 3 


Yes 


10" 


.2 


0.16 


0.075 


40.0* 


0.0068 


238.6 



Table 2. Results. The collision speeds \>n and V22 are averaged over orbits 20-30, while the column density ratio and 
boulder radius growth rate 02 (given here in Stokes number per orbit) are averaged over the last 10 orbits in the simulation. The 
recycling time-scale f recyc is calculated from Eqs. ( T37T i and ( BIT ). Measurements marked with * had not yet saturated at the end of 
the simulation. 
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regions must produce enough dust not only to feed its own 
zone, but also to fill up the regions that contain no boulders, 
as any gradients in the dust density will be quickly [instan- 
taneously actually, in the simplified model presented in Eqs. 
d30ll-(l3lTll evened out by diffusion. All together clumping leads 
maximally to a 25% increase in radius growth, but at the cost of 
reducing the total amount of boulders proportionally to the de- 
gree of clumping. We investigate the role of turbulent diffusion 
further in ^5] 



4.2. Recycling time-scale 

Even though a balance between sweep-up and collisional frag- 
mentation arises, so that the number density of fragments n\ 
stays approximately constant in time, there is a significant flux 
of dust grains through the boulder component. The evolution 
equation for «i [Eq. (Q3] consists of two terms that balance 
out in the equilibrium. One can rewrite the evolution equation 
in terms of the time-scale for sweep-up f swee p and the time- 
scale for replenishment of dust fragments by boulder collisions 

^replenish' 



1 dpi 1_ _J_ 

P 1 & t sweep ^replenish 



(34) 



In equilibrium the sweep-up and replenishment time-scales are 
equal. The sweep-up time-scale is 



1 



nzO~\2V\2 



(35) 



Inserting the equilibrium solution p\lp% = Avzzlvvz from 
Eq. Q gives 



1 + 4v 22 /v 12 a 2 p. 

Vl2 3pi +2 



(36) 



when assuming spherical grains with material density p.. One 
can simplify the equation further by inserting the Epstein 
regime friction time from Eq. (fT2l . yielding 



1 +4v 22 /vi 2 /pi +2 \ /Vn^ 1 



Tf 



(37) 



This is the time-scale upon which the boulders would empty the 
grain component in the absence of collisional fragmentation. 
In coagulation-fragmentation equilibrium collisional fragmen- 
tation produces dust grains at the same rate as they are swept 
up. But the sweep-up time-scale can also be associated with a 
characteristic recycling time-scale. On the average dust frag- 
ments spend the time f swee p in the grain component, before they 
are incorporated into a boulder. We have calculated the recy- 
cling time-scale for runs A and B, based on Eq. ( |37| |. in Table 
[2] The grains have a recycling time-scale of around three orbits 
in the two runs where dust fragments can not leave the boulder 
layer. Thus the grains have been through approximately 100 
agglomeration-destruction cycles during the course of runs A 
and B. We return to the recycling time-scale in models where 
dust fragments are allowed to diffuse out of the mid-plane layer 
in ED 



5. Diffusion 

In this section we generalise the analytical model of §|2] to 1- 
D. Adding a z-direction to the problem and exposing the dust 
particles to turbulent diffusion and vertical gravity yields the 
following equation system forpi = m\n\ and p2 = mini: 

dpi 
dt 



1 



1 



-pip 2 cr 12 Vi2 + p 2 cr 22 V 22 

m 2 m 2 



d(wipi) 
dz 



oz 



d(pi/p s ) 



dpi 
dt 



+PlP2C"i2Vi2- 



mi 



■ P 2 0"22V22- 



mi 



d(wipi) d 

H L>2 — 

dz dz 



dz 



d(pi/p s ) 



dz 



(38) 



(39) 



Here D\ and Di is the turbulent diffusion coefficient of grains 
and boulders, respectively, and w\ and wi are the vertical ve- 
locities, assumed to be in equilibrium between gravity and drag 
force with 



z 



Z ■ 



(40) 



This approximate expression recovers the terminal velocity of 
the grains in the small friction time regime, vv, = -TjQ^z, and 
the settling time of o scillating particles in the lar ge friction time 
regime, f sett = 1/r, (lYoudin & Lith wield, 120071) . The diffusion 
coe fficient D, depends on the friction time of the par ticles as 
(see lCarballido et allE006llYoudin & Lithwicki l2007h 

Do 



Di = 



(41) 



where Do is the diffusion coefficient of a passive scalar. We 
assume that the turbulent mixing is independent of the height 
over the mid-plane, i.e. that Do is a constant. Ignoring the frag- 
mentation and sweep-up terms of Eqs. (l38l and d39l allows for 
a simple equilibrium between diffusion and sedimentation, 



Pi(z) = 



■ eX p[-z 2 /(2Hf)], 



where the scale-height H, obeys the relation 
1 1 



H 2 + Do 



(42) 



(43) 



and 2Tj is the total column density of solids of type i. Here the 
denominators of Eqs. ( |40b and d4T1 ) have cancelled, making the 
short friction time scale height expression (Eq.[43l valid for all 
particle sizes. However, the equilibrium solution in Eqs. ( l42l - 
(l43l is formally only valid when the friction time r, is assumed 
constant with height over the mid-plane, an assumption that 
breaks down whe n considering the dust distribu tion over sev- 
eral scale heights (iDullemond & DominikL 120041) . 

Integrating Eqs. d38l and d39l over the entire z-space yields 
dynamical equations for the column densities Z\ and Si in- 
stead, 



dt ~ J_ a 

dx 2 r°°l 

dt ~ J A 



1 2 1 , 

-pip 2 cr 12 v 12 +p 2 cr 22 v 22 dz, 



mi 



mi 



1 2 1, 

+P1P2C12V12 p 2 cr 22 V 22 dz , 

\ mi mij 



(44) 
(45) 
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Fig. 5. Time series of the concentration of small dust grains for a run with a = 0.01 where the boulders lie in a thin layer around 
the mid-plane of the disc (run F). Fragments form when boulders collide in the mid-plane, but the produced dust grains escape 
from the mid-plane due to turbulent diffusion, spreading eventually evenly over the box. 



where the derivative terms have vanished because the column 
density of solids can not be changed by sedimentation and ver- 
tical diffusion. Inserting the sedimentation-diffusion equilib- 
rium solution from Eq. (l42l into Eqs. (144-b and d45T > and search- 
ing for E\ I Ei that gives E\ — E2 — yields 



E\ _ 0~22 V22 

£2 cr 12 V12 



\+H 2 jm 



V2 



(46) 



as an extension to Eq. (|5). If T! = T2, then H\ = H2 accord- 
ing to Eq. (l43l . and E\/E2 = (ct"22/o"12)(V22/vi2), completely 
equivalent to the 0-D case. Combining Eqs. d42b . d43l and d46b 
and inserting p\{z = 0) in Eq. (0]i yields the radius growth of 
boulders in the mid-plane as 



a.2 = 



E l+2 /(^H0 



V22 



(47) 



p. [2/(1 + #?/flf)] l ' 2 + 4V22/V12 ' 

the 1-D generalisation of Eq. ([6j. The growth rate of the boul- 
ders is more or less inversely proportional to Hi (note that the 



H\jH\ term in the denominator of Eq. l47l has only little in- 
fluence on the growth rate for H\ » #2). Thus going from 
Hi = H2 to, say, H\ = IO//2, a reasonable value for small frag- 
ments, decreases the growth rate by coagulation-fragmentation 
by a dramatic order of magnitude. 

The validity of the Gaussian solution (Eq. 11421 ) can be 
quantified by comparing the relevant time-scales of Eqs. ( |38| > 
and ( f39b - the time-scale for fragments to diffuse out of the 
mid-plane layer, r^iff = H^/D t , and the collisional time-scale, 
f coii = m 2/(P2C"22i'22)- The Gaussian solution is valid when 
fdiff •« fcoii, giving 



(48) 



where we used T2 = («2P.)/( c sPg) to translate the sedimen- 
tation time-scale into an Epstein friction time. Using fur- 
ther the sedimentation-diffusion equilibrium expression for 
the mid-plane density p2/p g = eo V£?KTf/£, where e is the 
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Fig. 6. Growth and fragmentation of boulders in models where the dust fragments are allowed to diffuse out of the boulder layer. 
The radius of the boulders increases much slower than in Fig.@] since fragments created in catastrophic collisions easily escape 
from the mid-plane layer and spread over the entire vertical extent of the disc. The ratio of fragments to boulders (panel 3) 
approximately matches the value expected from analytical derivations (panel 4) based on the different equilibrium scale height of 
boulders and fragments. Due to the slow growth boulders will have drifted into the inner disc to evaporate there before the meter 
barrier would eventually be crossed after a few thousand orbits. 



global solids-to-gas r atio of the boulder component (see e.g. 
Joha nsen et al. . 120061) . yields 



■1/3 



V2l\ 



2/3 



(49) 



with fully independent parameters (both the collision speed voi 
and the diffusion coefficient 5 of course depend on the strength 
of the turbulence, so these two must be chosen consistently). 
For 5 = 1CT 3 , eo = 0.01 and i'22/cs = 0.02 the criterion for 
the validity of the Gaussian solution is QkJi » 0.03, in ac- 
cordance with our modelling of component 2 as boulders with 

fi K Tf > 1. 

The simple model presented in 34. 1 I for the effect of clump- 
ing on coagulation-fragmentation also gives a new perspective 
on the effect of sedimentation [Eqs. (|46] i and ffll\\. Decreasing 
the scale height of boulders from H2 = Hi to H2 «: H\ leads 
to a steep rise in E1IZ2 [Eq. (|46|i1, but only to a mild increase 
in radius growth [Eq. (|4"7Ti1. Thus any non-homogeneity of the 
boulder layer (be it due to concentrations in transient gas high 



pressure or due to sedimentation) has little effect on the growth 
rate of the boulders, but may reduce the bulk density of boul- 
ders drastically. Taking instead the scale height of fragments 
Hi and increasing it from Hi = H2 to Hi » H2, a transi- 
tion which is similar to the one occurring from models A-B to 
models C-F, leads to a sharp decrease in both growth rate and 
column density of the boulder component. 



6. Results including leaking 

Having found in the preceding section an approximate analyt- 
ical solution for the coagulation-fragmentation equilibrium in 
the case where dust fragments are free to leave the sedimentary 
boulder layer, we return now to the results of numerical sim- 
ulations. In runs C-F we set the solids-to-gas ratio of boulders 
and fragments to the more canonical value 0.01 and expose the 
boulders to vertical gravity. We give the boulders time to set- 
tle to the mid-plane from t = 10r m b to t — 20r ol -b, so that an 
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equilibrium sedimentary layer has formed when sweep-up and 
collisional fragmentation are turned on at t — 20r or [,. 

In Figs. [3] and [6] we show the time evolution of fragments 
and boulders in models where the boulders have sedimented 
out of the gas to establish a thin layer around the mid-plane 
of the box. Fragments form in catastrophic collision between 
boulders in the mid-plane, but are quickly carried away from 
the mid-plane layer by the turbulent gas and are eventually 
well-mixed throughout the box (Fig. 0. Fig. [6] shows, for two 
different numerical resolutions, that the boulders grow a fac- 
tor ten slower than in the case where dust fragments were not 
allowed to leave the boulder layer (Fig. 0). Also the mass of 
fragments is huge because the catastrophic collisions have to 
keep up with the turbulent transport of fragments away from 
the mid-plane. The fourth panel of Fig. [6] shows the analytical 
expectation value for M\ IM%, following Eq. (|46*| |. by setting the 
expected scale height of fragments, H\ , equal to the gas scale 
height and using the collision speeds from the second panel. 
There is a bit more mass in fragments than expected, which 
may, as in the case where dust fragments were not allowed to 
leave the boulder layer shown in Fig. |4] be explained by the 
clumpy structure of the boulder layer (see §14.11 ). 

The turbulent transport reduces the column density of boul- 
ders to approximately 1-2% of the total column density of the 
solids. The growth rate of the Stokes number of the boulders is 
around 0.003 . . . 0.007 per orbit, with the higher values appear- 
ing in the strongly turbulent a = 10~ 2 runs. Under all circum- 
stances this is way too low to compete with radial drift which 
occurs on a time-scale of a few ten orbits. Even if the radial drift 
time-scale is ignored it would take around 1000 orbits to grow 
to St = 10. The extremely small column density of the boulders 
corresponds to around 1 boulder per (10 km) 2 . It is somewhat 
surprising that so few boulders can popula te the entire vertica l 
extent of the disc with small dust grains. Eisner et al.l (|2006) 



recently modelled submicron-sized dust grains in the inner part 
of the transition disc TW Hya and found that the lifetimes of 
these grains against radiation pressure is much shorter than the 
age of the system. One can speculate that the source of these 
small grains could be boulders (drifting through the inner disc 
or permanently situated there) creating observable amounts of 
fragments as they collide. 

The minimum and maximum bulk density of fragments in 
the gas is shown in Fig. [7] After collisional fragmentation is 
turned on at t — 20r ol -b, there is a sharp increase in the maxi- 
mum concentration of fragments, but turbulent mixing eventu- 
ally leads to a state where the minimum and maximum concen- 
trations are within 3% of each other for a = 10~ 2 and within 
6% of each other for a = 10~ 3 . It may be surprising that the dif- 
ference is so little, but it is again because the collisional time- 
scale is much longer than the diffusion time-scale. That also ex- 
plains why the a = 10~ 2 has less difference between minimum 
and maximum concentrations than a = 10~ 3 - the stronger 
mixing in the highly turbulent case evens out concentration dif- 
ferences more effectively. We note that although the average 
solids-to-gas ratio in the disc is assumed to be 6q = 0.02, the 
fragments reach an average bulk density of 6q = 0.029. This is 
an artificial effect of the limited vertical extent of the box: boul- 
ders have sedimented out of the regions outside of the box, but 
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Fig. 7. The minimum and maximum bulk density of fragments, 
pi, in the box as a function of time for runs E and F. Collisional 
fragmentation is turned on at t — 20r or b, followed by an initial 
peak in the maximum density of fragments. Turbulent mixing 
nevertheless leads to an almost homogeneous state with only 
very little difference between minimum and maximum con- 
centrations of fragments. The case a = 10 -2 has less differ- 
ence between minimum and maximum concentrations than for 
or = 10 -3 due to stronger mixing in the first case. 



fragments are not allowed to escape back out the box, hence 
the average density of fragments is kept artificially high. 

6.1. Recycling time-scale revisited 

We introduced the recycling time-scale in 94.2I With diffusion 
of dust fragments from the boulder layer suppressed, the dust 
grains would be free for only around three orbits before being 
incorporated into a boulder. We now derive a similar expres- 
sion for the recycling time-scale, valid in the case where dust 
fragments are allowed to diffuse out of the thin mid-plane layer 
where they originate in collisions between boulders. 

Assuming that the density dependence on height over the 
mid-plane is Gaussian for both species (Eq. l4"2ll ) we can eval- 
uate the integral in Eq. ( PPfl i analytically and get the recycling 
time-scale 



m 2 



^recyc 



^2 



C12V12 



(50) 



Notice that the recycling time-scale for H\ = Hi is V2 times 
longer than in the case where dust fragments and boulders were 
not allowed to separate vertically, because, although the densi- 
ties in the very mid-plane are exactly as in the 0-D case, the 
vertically averaged collision rate is smaller, which leads to a 
somewhat longer recycling time-scale. 
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Inserting the coagulation-fragmentation equilibrium from 
Eq. (|46| i yields the recycling time-scale in units of the friction 
time as 



1 + 4(v 22 /v 12 ) ( ^1 + H\IH\I V5) 



Tf 



Vl2 



Cs 



Pg 



(51) 



completely equivalent to Eq. ( f3Tb . Using S\+2 — Q.Q2£ S = 
0.02 V27T in units where the mid-plane gas density is unity, to- 
gether with H\/H g = 1 and the diffusion-sedimentation equi- 
librium Hx/Hg = y/5/(QKTf), we have calculated the recycling 
time-scale based on Eq. ( f5Tb in Table [2] Allowing dust frag- 
ments to leak out of the mid-plane layer leads to a dramatic in- 
crease by two orders of magnitude in the recycling time-scale, 
which is now 200-300 orbits for the a = 10~ 2 models and ap- 
proaching a thousand orbits for the a = 10~ 3 runs. Not only 
does the vertical diffusion decrease the amount of boulders in 
the mid-plane by approximately an order of magnitude, leading 
to a much longer sweep-up time-scale, but the total amount of 
dust grains is also higher, so that it takes much longer time for 
the average dust grain to encounter a boulder. 

7. Summary and discussion 

We have proposed a simple two-component model for the 
growth and collisional destruction of boulders in protoplane- 
tary discs. Fragments produced in catastrophic collisions be- 
tween boulders are swept up by other boulders, leading to a 
continuous growth towards larger bodies. An analytical equi- 
librium solution to the dynamical equations predicts that the 
boulder radius can grow as quickly as a few mm per year. 
The promise of this method to provide growth rates that can 
compete with the radial drift is nevertheless compromised by 
the inclusion of turbulent transport of fragments out of the 
boulder layer. The sedimentary mid-plane layer loses the pro- 
duced fragments to the boulder-free parts of the disc, eventually 
grinding down the boulder component to insignificant masses 
and reducing the growth rate drastically. Thus what initially 
seemed to have the potential to provide an efficient growth 
phase of boulders turned into something like a worst-case sce- 
nario with the inclusion of turbulent transport. 

One can think of ways by which boulders may still pene- 
trate the meter barrier by sweep-up of small fragments: 

- Fragments are larger 

- Boulder collisions do not lead to destruction 

- Protoplanetary discs are less turbulent than assumed 

- Turbulence is confined to the mid-plane 

- Radial drift is weaker in nature than thought 

If fragments are large enough that they do not couple in- 
stantaneously (compared to an orbital time-scale) to the gas, 
then the turbulent transport away from the sedimentary mid- 
plane layer is slowed down. Having such large fragments 
would in principle put coagulation-fragmentation growth back 



on track. BDH presented models with a much more advanced 
collisional fragmentation model, where the results of catas- 
trophic collisions are distributed with a power law, but this 
still did not lead to a break-through of the meter barrier. If, 
on the other hand, boulders do not fragment at all, but merely 
bounce off each other, then the whole concept of coagulation- 
fragmentation growth breaks down, and the stage is left to 
coagulation of equa l-sized boulders and/o r self-gravity in the 
boulder component jjohansen et all 120071) . 

Turbulence plays an interesting double role in the 
coagulation-fragmentation process. The relative speed between 
the boulders, which leads to their destruction and the contin- 
uous replenishment of the grains, is induced by the marginal 
coupling of the boulders to the turbulent gas motion. On the 
other hand turbulent diffusion drains the mid-plane layer of 
fragments and reduces the growth rate of the boulders. This 
is in some opposition to the role that is normally attributed 
to turbulence in the coagulation process. Here a dense mid- 
plane layer and reduced collisional fragmentation are desired, 
and turbulence counteracts both. Coagulation-fragmentation 
growth, on the other hand, benefits directly from stronger tur- 
bulence and higher collision speeds (although the sweep-up ef- 
ficiency can be put into question when the relative speed be- 
tween dust grains and boulders incr eases beyond a few t en me - 
ters per second, see discussions in IWurm et all 1200 U 120051) . 
while it is indifferent to sedimentation and clumping of the 
boulder layer. 

The presented simulations all have a simplified space- 
filling turbulence of magnetic origin. If turbulence would in- 
st ead be confined t o the mid-plane of the disc, as assumed 
in Weidenschillind ( 1997 ). then the transport of grains away 
from the mid-plane can be reduced, while at the same time 
the collision speeds can be kept high. The Kelvin-Helmholtz 
and streaming instabilities associated with the sedimentation 
of solids are nevertheless not Keplerian shear instab ilities 
( lYoudin & Goodmarl 12005c lYoudin & Johansenl 120071) and 
thus cannot explain the observed accretion rates of young 
stars. If, on the other hand, there is a region around the mid- 
plane where the magnetic field does not coupl e well enough to 
the gas to have m agnetorotational instability (IGammie , 1 1996c 



Oishi et al. , 2007), then one could have lower turbulence in 



the mid-plane and much less loss of fragments to the boulder- 
free parts of the disc, but as discussed above, a net decrease 
of turbulence in the mid-plane has a negative impact on the 



coagulation-fragmentation growth. Dust gra ins are a ma 



alyst for recombination of ionised species ( ISano et al 



orcat- 
20001) . 



This leads to an interesting coupling between dust and turbu- 
lence whereby dusty regions would have weaker turbulent mo- 
tion, and thus less diffusion, than dust-free regions. We plan to 
include such effects in a future model. 

Collisional fragmentation is actually not the real problem 
for planetesimal formation, given that fragments are readily in- 
corporated in the few lucky boulders that avoid catastrophic 
collisions. The problem arises because the time-scale of growth 
by coagulation-fragmentation is so long that all material will 
have been flushed through the disc long before being able to 
grow big enough to decouple from the gas. Radial drift in the 
minimum mass solar nebula reaches 5% of the local sound 
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speed ( Weidenschillind, 1977 ). yielding a drift time-scale of a 
few ten orbits, much too short for any significant radius growth. 

There is a mounting observational evidence that radial drift 
may not be as big in actual discs as predicted from theoreti- 
cal arguments. iRettig et al. I d2006l) measured the abundance oi" 
gas and dust in discs of millions of years in age and were able 
to explain the inclination dependence of the dust-to-gas ratio 
from a combination of grain growth and sedimentation, still 
within the framework of the minimum mass solar nebula with 
a solids-to-gas ratio of 0.01. In the models of drift and coag- 
ulation presented in BDH, on the other hand, the outer disc is 
cleared of dust in a few hundred thousand years. The observed 
presen ce of cm-sized solid s in the oute r parts of protoplanetar 



discs dWilner et all 120001 iTesfi et all 120031: iRodmann et" 



ary 
al., 



20061: Lommen et all l2007l)~ is at best marginally consistent 



with theoretical lif e-times of such grains due to radial drift 
dBrauer et alll2007l) . Maybe the most promising way to stop ra- 
dial d rift is to have radial pressure bumps in the disc ( Whipple, 
1972ft . These bumps can arise from first principles in 3-D 
simulations of the dynamics of pr otoplanetary discs , e.g. spi- 



ral arms of self-gravitating discs ( IRice et all 120041) or long 



lived high pressure regions in magnetoro t ationa l turbulence 



dFromang & Nelsonl l2005t iJohansen et al. . 120061) . Any pres 



sure bump must compete with the global pressure gradient and 
produce a net zero or positive gradient. This requirement is nev- 
ertheless not necessarily very difficult to fulfil. The gas pressure 
typically falls 10% over one scale height in the radial direc- 
tion, a value with which even subsonic turbulence can easily 
compete. The challenge is to have long-lived pressure enhance- 
ments, an issue which is still debated for the cas e of magne- 
torotational turbulence dFromang & Nelsonl 120061) . Radial drift 
may also be significantly reduced in a dense boulder-dominated 
mid-plane layer where the gaseous he adwind is reduced as th e 



gas is dragged along with the boulders ( Nakagawa et al. , 19861) 



The local absence of radial drift would not only solve the 
time-scale problem, but also reduce differential radial drift, 
which was shown by BDH to lead to destruction of the boulders 
due to collisions at speeds between 10 and 50 m/s with slightly 
smaller boulders. The process of crossing from 1 m to 10 m 
would still be very inefficient, since the turbulent transport 
would cause 99% of the solid mass to be present in small frag- 
ments, but given that planet formation in our solar system took 
place over millions of years (see review by iTrieloff & Palme , 
20061) . inefficient planetesimal formation may actually be de- 
sired to comply with meteoritic evidence. Once a few extraor- 
dinarily lucky bodies would cross the meter barrier, they could 
grow big enough to even sweep up the boulders still lying 
around in the mid-plane and thus continue to grow towards 
full-fledged protoplanets and later gas giant cores and terres- 
trial planets. 
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Appendix A: Stability analysis 

In this appendix we consider the stability of the equilibrium 
solution to Eq. ([T). We find that the amplitude of any (arbi- 
trarily large) perturbation to the equilibrium will decrease ex- 
ponentially with time, so that Eq. d5j is a (both linearly and 
non-linearly) stable solution to the coagulation-fragmentation 
problem. 



We start by rewriting Eq. (Q3 slightly by multiplying with 
mi and inserting pi + 2 — p\ + pi (mass conservation) to avoid 
any reference to pi- The resulting dynamical equation is 



dpi P\+2~P\ (Pi+2-Pi) 2 

-Z- = "Pi 0- n V\2 + cr 22 V 2 2 • 

at rri2 m.2 



(A.l) 



We linearise this equation by writing the density of fragments 
Pi as 



Pi = Pi +P\(t), 

where p l is the equilibrium solution to Eq. dA.l 



Pi =Pl+2 



C 22^22 



c 12V12 + cr 22^22 



(A.2) 



(A.3) 



and p'j is an infinitesimal perturbation to this equilibrium den- 
sity. We ignored the other equilibrium solution, p l = pi+2, 
because that state is not accompanied by any radius increase 
of the boulders (see discussion below). The expression in 
Eq. dA.3b is completely similar to Eq. (0, but we avoided 
approximating cr 2 2 ~ 4o"i2 for generality reasons. Inserting 
Eq. (lA.2b into Eq. (IA.U yields the dynamical equation 



d P\ ^ ,, Pi+2-Pi -P\ 

— = ~{p l +p l ) <T 12 V12 

at ni2 

(Pl+2-Pl -P'i) 2 
+ - 

rri2 



-cr 22 V22 



(A.4) 



for the density perturbation p' x . Using the fact that p x in itself 
satisfies Eq. ( IA.1I ) and ignoring terms of second order in the 
perturbed density yields the linearised dynamical equation 



dp[ 
dt 



Q~\2V\2 

m 2 



(- 2 PlPl +p'iPl+2) 



+ ^22V22 ( _ 2pi+2p , + 2p lP \). (A.5) 

m.2 



Inserting now the equilibrium solution ~p x from Eq. ( IA.3b gives, 
after some trivial algebraic manipulation, the final equation for 
the evolution of the perturbed density as 



dt 



Pl+2 



-O-12V12P1 



(A.6) 



Thus any small variation from the equilibrium expression will 
decay exponentially [following p' x (t) oc exp(-wf)] at the rate 



Pl+20" 12V12 



(A.7) 



The corresponding decay time W2/(pi+2t r i2Vi2) is, not surpris- 
ingly, similar to the time-scale of the sweep-up process, but 
with the the density of the fragments replaced by the total den- 
sity pi +2 . 

One may suspect that the equilibrium solution is also non- 
linearly stable, due to the simple nature of the Eq. ( lA.ll i. 
Rewriting Eq. ( IA. lb in terms of the normalised density pi = 
P1/P1+2, normalised time t = f/[m2/(pi+20"i2Vi2)] and collision 
parameter £ = cr22V22/(c r i2Vi2) yields the evolution equation 



^=(l+£)p 2 +(-l-2£)p 1+ £. 



(A.8) 
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Fig.A.l. The dependence of the time derivative of the nor- 
malised density of fragments, dpi/dt, on the density, p\ = 
pi/pi+2, itself, following Eq. dA.8b . All states tend towards the 
equilibrium state at p\ = £/(l + £). Here £, = (TiivziKcr^vn) is 
a parameter that depends on the collision speeds and on the 
cross section of the boulders. The equilibrium at pi = 1 is 
clearly unstable. 



The second order polynomial dpi/dt = has the solutions 



n = l, 



ri = 



(A.9) 
(A.10) 



In the first solution all the solid mass is bound in frag- 
ments - the total absence of boulders for this case means that 
coagulation-fragmentation growth is disabled. The second so- 
lution is the same as in Eq. dA.3l ). Since we know that the poly- 
nomial in Eq. d A.8I > opens upwards and that there is only one 
crossing of zero in the interval [0, 1[, then all states of p\ in 
this interval must approach the equilibrium solution given in 
Eq. (IA. 10b . See Fig. IA.1I for an illustration. The same argu- 
ments also show that p\ — 1 [Eq. dA.9H is an unstable solution, 
since even a vanishingly low number density of boulders will 
send the state towards ri instead. One must nevertheless still 
require some minimum number of boulders in the system for 
the continuity description of their number density to hold. 



using spline interpolation (see Youdin & JohansenT. 2007 ). The 
Epstein drag law, with the friction time tf given by 



(B.l) 



is valid as long as the particle radius «2 < (9/4)/l, where A is 
the mean free path of the gas molecules. At «2 > (9/4)/i one 

(St) 

needs to use instead the Stokes friction time, given by Qkt { = 
^KTf p) x(4/9)(a 2 /^).The mean free path of the gas molecules 
can be calculated from 



A = 



Pg°"mol ^gCmol 



I 1 



(B.2) 



where [i = 3.9 x 10 g is the mean molecular weight and 
2.0 x 10~ 15 cm 2 is the cross section of molecular h 



c mo i 



drogen dChapman & Cowling , Il970fc iNakagawa et all Il98 
In units of the gas scale height H the mean free path can be 
expressed as 



,1 
H 



4.9 x lO^gcirT 



(B.3) 



Assuming a minimum mass solar nebula model at r — 5 AU 
we have £ g = 150 gem 2 and therefore a mean free path of 
A/H = 3.3 x 10~ n . The transition from Epstein to Stokes 
regime thus occurs at azIH = 7.3 x 10~". We start our boul- 
ders with radius ao/H = 10~ n (with U^Tf = 1), so the Epstein 
regime is valid up to £?KTf ~ 7.3, which is already at the other 
side of the meter barrier. Since the focus of this paper is the 
crossing of the meter barrier, we shall for simplicity model drag 
force in the Epstein regime throughout and ignore the transition 
to the Stokes regime. The coagulation-fragmentation equilib- 
rium anyway has no dependence on the assumed amount of 
gas in the disc - the ratio of dust fragments to boulders de- 
pends only on the collision speeds [Eq. (0], and while the ra- 
dius growth [Eq. ©J does scale with pi+2, and thus with p g if 
the solids-to-gas ratio is unchanged, the evolution of the Stokes 
number, St oc ailpg, is independent of the gas density. 

Another complication with modelling the Stokes regime 
is that the gas flow in the vicinity of a boulder would lead 
small dust grains around the boulder rather than onto its sur- 
face (ISekiva & Takedal 120031) . Dust grains may however still 
be able to penetrate to the b oulder in case the bo ulder is porous 



and has gas flow through it dWurm et al.L 120041) . We shall nev 



ertheless limit ourselves to the Epstein regime in this paper and 
leave the treatment of the interaction of grains and boulders in 
the Stokes regime to a future publication. 



Appendix B: Drag force 

In this appendix we describe the implementation of drag forces 
in our simulations. We let gas exert drag on the boulders fol- 
lowing an Epstein drag law that is linear in the velocity differ- 
ence between particles and gas. The gas velocity at the posi- 
tion of a particle is interpolated from the 27 nearest grid points 



